import Data.List
import Control.Monad
import Data.Maybe
import Data.IntMap (IntMap)
import qualified Data.IntMap as IntMap

ourPrimes = takeWhile (<10000) primes

--pairs :: IntMap [Int]
--pairs = IntMap.fromList [ (prime,makePairs prime) | prime <- ourPrimes ]
pairsList :: [(Int, Maybe [Int])]
pairsList = zip ourPrimes (map makePairs' ourPrimes)

makePairs' :: Int -> Maybe [Int]
makePairs' prime
	| pairList == [] = Nothing
	| otherwise = Just pairList
	where
	pairList = [p | p <- dropWhile (<prime) ourPrimes, checkConcat prime p]

makePairs :: Int -> Maybe [Int]
makePairs x = fromJust $ lookup x pairsList
	
checkConcat :: Int -> Int -> Bool
checkConcat p1 p2 = (isPrime $ read ((show p2) ++ (show p1)) ) &&
						(isPrime $ read ((show p1) ++ (show p2)))

						
nextValues :: ([Int],Maybe [Int]) -> [([Int],Maybe [Int])]
nextValues ( startList , Nothing ) = []
nextValues ( startList , Just [] ) = []
nextValues ( startList , Just nextList ) = 
	map getNextList nextList
	where
		getNextList x = (startList ++ [x],liftM2 intersect (Just nextList) (makePairs x))
			
newLists = [ [a,b,c,d] | a <- primes, b <- takeWhile (<a) primes, c <- takeWhile (<b) primes, d <- takeWhile (<c) primes]

testNewList testList = and [ isPrime $ combine a b | a<-testList, b<-testList, a /= b ]
	where
		combine a b = read $ show a ++ show b
		
winningList = head $ dropWhile (not.testNewList) newLists

euler_060 = (sum.fst.head) $ nextValues ([],Just ourPrimes) >>= nextValues >>= nextValues >>=nextValues >>=nextValues

-- "Library" Functions below

primes :: Integral a => [a]
primes = 2 : filter ( (==1) . length . primeFactors ) [3,5..]

primeFactors :: Integral a => a -> [a]
primeFactors n = factor n primes
	where
	factor n (p:ps)
		| p*p > n		 = [n]
		| n `mod` p == 0 = p : factor ( n `div` p )	(p:ps)
		| otherwise		 = factor n ps

isPrime x
	| x==3=True
	| otherwise = and [ millerRabinPrimality x n | n<-[2,3] ]

-- (eq. to) find2km (2^k * n) = (k,n)
find2km :: Integral a => a -> (a,a)
find2km n = f 0 n
    where 
        f k m
            | r == 1 = (k,m)
            | otherwise = f (k+1) q
            where (q,r) = quotRem m 2        
 
-- n is the number to test; a is the (presumably randomly chosen) witness
millerRabinPrimality :: Integer -> Integer -> Bool
millerRabinPrimality n a
    | a <= 1 || a >= n-1 = 
        error $ "millerRabinPrimality: a out of range (" 
              ++ show a ++ " for "++ show n ++ ")" 
    | n < 2 = False
    | even n = False
    | b0 == 1 || b0 == n' = True
    | otherwise = iter (tail b)
    where
        n' = n-1
        (k,m) = find2km n'
        b0 = powMod n a m
        b = take (fromIntegral k) $ iterate (squareMod n) b0
        iter [] = False
        iter (x:xs)
            | x == 1 = False
            | x == n' = True
            | otherwise = iter xs
 
-- (eq. to) pow' (*) (^2) n k = n^k
pow' :: (Num a, Integral b) => (a->a->a) -> (a->a) -> a -> b -> a
pow' _ _ _ 0 = 1
pow' mul sq x' n' = f x' n' 1
    where 
        f x n y
            | n == 1 = x `mul` y
            | r == 0 = f x2 q y
            | otherwise = f x2 q (x `mul` y)
            where
                (q,r) = quotRem n 2
                x2 = sq x
 
mulMod :: Integral a => a -> a -> a -> a
mulMod a b c = (b * c) `mod` a
squareMod :: Integral a => a -> a -> a
squareMod a b = (b * b) `rem` a
 
-- (eq. to) powMod m n k = n^k `mod` m
powMod :: Integral a => a -> a -> a -> a
powMod m = pow' (mulMod m) (squareMod m)
